rm(list = ls())

#setwd("~/Dropbox/Energiewende_Project/All_Replication_Files")
setwd("C:/Users/chris/Dropbox/Energiewende_Project/All_Replication_Files")

#Loading in Packages ----------------------------------------------------------------------------------

library(dplyr)
library(sandwich)
library(lmtest)
library(margins)
library(ggplot2)
library(scales)
library(gridExtra)
library(readr)
library(stargazer)
library(reshape2) 
library(stringr)
library(cregg)
library(kableExtra)
library(knitr)
library(modelsummary)
options(knitr.table.format = "latex")


####################################################################################################################################
############################################Policy Plan Experiment###########################################################
####################################################################################################################################

#Data Prep -----------------------------------------------------------------------------------------------------
data <- read_csv("data.coal.cleaned.final.csv", col_types = cols(.default = "?"))
data_wtp <- read_csv("data.wtp.final.csv", col_types = cols(.default = "?"))
weights <- read_csv("weights_data.csv")

data_wtp <- merge(data_wtp, weights, by=c("tic"), all.x=T)
data <- merge(data, weights, by=c("tic"), all.x=T)

#data prep
data_wtp$cost <- relevel(as.factor(data_wtp$cost), 
                         "150 EUR per Year")

data_wtp$effectiveness <- relevel(as.factor(data_wtp$effectiveness),
                                  "25% CO2 Reduction by 2030")

data_wtp$compensation <- relevel(as.factor(data_wtp$compensation),
                                 "No Change")

data_wtp$competition <- relevel(as.factor(data_wtp$competition),
                                "No Change")

data_control <- subset(data_wtp, data_wtp$control==1)

data_wtp$vignettetreatment <- as.factor(data_wtp$willingnesstopay)


########helper functions########
#regression functions
wtp_baseline_reg <- function(dataset, subgroup){
  dats <- dataset[dataset[,subgroup] == 1,]
  
  reg <- lm(choice_plan ~ cost + effectiveness + compensation + competition, data=dats, weights = dats$weights)
  
  vcreg <- vcovCL(reg, cluster=dats$ResponseId)
  reg_rob <- coeftest(reg, vcov=vcreg)
  return(reg_rob)
}

wtp_baseline_reg_contr <- function(dataset, subgroup){
  dats <- dataset[dataset[,subgroup] == 1,]
  
  reg <- lm(choice_plan ~ cost + effectiveness + compensation + competition +
              region + incomequintiles + gender + age, data=dats, weights = dats$weights)
  
  vcreg <- vcovCL(reg, cluster=dats$ResponseId)
  reg_rob <- coeftest(reg, vcov=vcreg)
  return(reg_rob)
}

wtp_interactions <- function(dataset, subgroup, interaction){
  
  dats <- dataset[dataset[,subgroup] == 1,]
  
  xvars <- paste0("cost + effectiveness + compensation + competition+ ",interaction,"+ compensation*",interaction ," + competition*",interaction)
  yvar <- "choice_plan"
  fmla.reg = as.formula(paste0(yvar, '~', xvars))
  
  regout <- lm(fmla.reg, data = dats, weights = dats$weights) 
  
  
  vcreg <- vcovCL(regout, cluster=dats$ResponseId)
  reg_rob <- coeftest(regout, vcov=vcreg)
  return(reg_rob)
}

wtp_interactions_contr <- function(dataset, subgroup, interaction){
  
  dats <- dataset[dataset[,subgroup] == 1,]
  
  xvars <- paste0("cost + effectiveness + compensation + competition+ ",interaction,"+ compensation*",interaction ," + competition*",interaction,"+region + incomequintiles + gender + age")
  yvar <- "choice_plan"
  fmla.reg = as.formula(paste0(yvar, '~', xvars))
  
  regout <- lm(fmla.reg, data = dats, weights = dats$weights) 
  
  
  vcreg <- vcovCL(regout, cluster=dats$ResponseId)
  reg_rob <- coeftest(regout, vcov=vcreg)
  return(reg_rob)
}

wtp_interactions_robust <- function(dataset, subgroup, interaction){
  
  dats <- dataset[dataset[,subgroup] == 1,]
  
  xvars <- paste0("cost + effectiveness + compensation + competition+ ",interaction,"+ compensation*",interaction ," + competition*",interaction,  "+ climatehumans*compensation + climatehumans * competition + left_right_4*compensation + left_right_4*competition + party_2021_2*compensation + party_2021_2*competition + climate_policyn*compensation + climate_policyn*competition")
  yvar <- "choice_plan"
  fmla.reg = as.formula(paste0(yvar, '~', xvars))
  
  regout <- lm(fmla.reg, data = dats, weights = dats$weights) 
  
  
  vcreg <- vcovCL(regout, cluster=dats$ResponseId)
  reg_rob <- coeftest(regout, vcov=vcreg)
  return(reg_rob)
}

#functions to draw plots

data.prep.wtp <- function(reg){
  #get rid of constant
  reg <- as.data.frame(reg[2:11,])
  reg$factor <- rownames(reg)
  
  #add the comparison categories:
  reg[11:14,] <- 0
  reg$factor[11:14] <-   c("competitionNo Change", "compensationNo Change", "effectiveness25 Percent CO2 Reduction by 2030", "cost150 EUR per Year")  
  reg$factor <- as.factor(reg$factor)
  
  #introduce factor levels for different variables:
  reg$factor <- factor(reg$factor, levels = c(levels(reg$factor),
                                              "Competitiveness Policy",
                                              "Compensation Policy",
                                              "Effectiveness",
                                              "Cost"))
  reg$factor <- factor(reg$factor,
                       levels = levels(reg$factor)[
                         c(17,13,14, #effectiveness
                           18,10,11,12,9, #cost
                           16,5,3,1,4,2, #compensation
                           15,8,6,7 #competition policy
                         )
                       ])
  reg$factor <- factor(reg$factor, levels = rev(levels(reg$factor)))
  reg$confintlower <- reg$Estimate - 1.96* reg$`Std. Error`
  reg$confintupper <- reg$Estimate + 1.96*reg$`Std. Error`
  return(reg)
}

data.prep.wtp.mm <- function(mms){
  #get rid of constant
  reg <- mms %>% select(estimate, lower, upper)
  reg$factor <- c("cost150", "cost1200", "cost300", "cost600",
                  "25red", "55red", "compnochange", "comphighlumpsum",
                  "comphighprog", "complowlump", "complowprog", 
                  "competcbt", "competptas", "competnochange" 
  )
  
  #add the comparison categories:
  
  reg$factor <- as.factor(reg$factor)
  
  #introduce factor levels for different variables:
  reg$factor <- factor(reg$factor, levels = c(levels(reg$factor),
                                              "Competitiveness Policy",
                                              "Compensation Policy",
                                              "Effectiveness",
                                              "Cost"))
  reg$factor <- factor(reg$factor,
                       levels = levels(reg$factor)[
                         c(17,1,2, #effectiveness
                           18,12,13,14,11, #cost
                           16,10,8,6,9,7, #compensation
                           15,4,3,5 #competition policy
                         )
                       ])
  reg$factor <- factor(reg$factor, levels = rev(levels(reg$factor)))
  return(reg)
}
plot_conjoint <- function(regresults, title, filename){
  
  plotdata <- data.prep.wtp(regresults)
  pl <- ggplot(data=plotdata)+
    geom_point(aes(x = Estimate, y = factor))+
    geom_segment(aes(x=confintlower, xend = confintupper,
                     y = factor, yend = factor))+
    geom_vline(xintercept = 0, lty = "dotted")+
    theme_bw()+
    scale_y_discrete(NULL, drop=FALSE,
                     labels = rev(c("Effectiveness",
                                    "25% CO2 Reduction by 2030",
                                    "55% CO2 Reduction by 2030",
                                    "Cost",
                                    "150EUR per Year",
                                    "300EUR per Year",
                                    "600EUR per Year",
                                    "1200EUR per Year",
                                    "Compensation Policy",
                                    "No Change",
                                    "Low Lump Sum Carbon Dividend",
                                    "High Lump Sum Carbon Dividend",
                                    "Low Progressive Carbon Dividend",
                                    "High Progressive Carbon Dividend",
                                    "Competitiveness Policy",
                                    "No Change",
                                    "Carbon Border Tax",
                                    "Environmental Clauses in PTAs")))
  
  pl <- pl + theme(axis.text.y = element_text(
    face = rev(c('bold', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain','plain',
                 'bold', 'plain', 'plain', 'plain'))
  ))
  
  pl <- pl + ggtitle(title)
  pl <- pl + theme(text = element_text(size=16),
                   axis.text.y = element_text(size=16),
                   plot.title = element_text(size=18, hjust = 0.5),
                   axis.title.x=element_text(size=16)
  )+xlab("Average Marginal Component Effects on Support for Policy Plan")
  
  
  print(pl)
  
  
  ggsave(filename = paste0("Graphs/weighted",filename ,".png"),
         plot = pl,
         width = 12, height = 8)
  return(pl)
  
}

splitsampleplots_wtp <- function(datasets, splitvar, titleadd, filename){
  dataset <- datasets
  reg_wtp_hs_c <- lm(choice_plan ~ cost + effectiveness + compensation + competition, data=dataset[dataset[,splitvar]==1,])
  
  vcreg_wtp_hs_c <- vcovCL(reg_wtp_hs_c, cluster=dataset[dataset[,splitvar]==1,]$ResponseId)
  reg_wtp_hs_cr  <- coeftest(reg_wtp_hs_c, vcov=vcreg_wtp_hs_c)
  
  reg_wtp_ls_c <- lm(choice_plan ~ cost + effectiveness + compensation + competition, data=dataset[dataset[,splitvar]==0,])
  
  vcreg_wtp_ls_c <- vcovCL(reg_wtp_ls_c, cluster=dataset[dataset[,splitvar]==0,]$ResponseId)
  reg_wtp_ls_cr  <- coeftest(reg_wtp_ls_c, vcov=vcreg_wtp_ls_c)
  
  plot.wtp_hs_cr <- data.prep.wtp(reg_wtp_hs_cr)
  plot.wtp_ls_cr <- data.prep.wtp(reg_wtp_ls_cr)
  
  plot.wtp_hs_cr$`Belief in State Intervention` = paste0("High, n=", sum(data[,splitvar]==1, na.rm=T))
  plot.wtp_ls_cr$`Belief in State Intervention` = paste0("Low, n=", sum( data[,splitvar]==0, na.rm=T))
  
  plot.wtp_c_si <- rbind(plot.wtp_hs_cr, plot.wtp_ls_cr)
  
  plgs <- ggplot(plot.wtp_c_si)+
    geom_point(aes(x = Estimate, y = factor, color=`Belief in State Intervention`, shape = `Belief in State Intervention`), position=position_dodge(width=.75))+
    geom_linerange(aes(xmin=confintlower, xmax = confintupper,
                       y = factor,  color=`Belief in State Intervention`), position=position_dodge(width=.75))+
    geom_vline(xintercept = 0, lty = "dotted")+
    theme_bw()+
    scale_y_discrete(NULL, drop=FALSE,
                     labels = rev(c("Effectiveness",
                                    "25% CO2 Reduction by 2030",
                                    "55% CO2 Reduction by 2030",
                                    "Cost",
                                    "150EUR per Year",
                                    "300EUR per Year",
                                    "600EUR per Year",
                                    "1200EUR per Year",
                                    "Compensation Policy",
                                    "No Change",
                                    "Low Lump Sum Carbon Dividend",
                                    "High Lump Sum Carbon Dividend",
                                    "Low Progressive Carbon Dividend",
                                    "High Progressive Carbon Dividend",
                                    "Competitiveness Policy",
                                    "No Change",
                                    "Carbon Border Tax",
                                    "Environmental Clauses in PTAs")))
  
  plgs <- plgs + theme(axis.text.y = element_text(
    face = rev(c('bold', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain','plain',
                 'bold', 'plain', 'plain', 'plain'))
  ))
  
  plgs <- plgs + ggtitle(paste0("Treatment Effects Policy Plan Conjoint,  \n by Belief in State Intervention, Weighted Regressions", titleadd))+
    theme(text = element_text(size=16),
          axis.text.y = element_text(size=16),
          
          plot.title = element_text(size=18, hjust = 0.5),
          legend.position = "bottom",
          axis.title.x=element_text(size=16),
          legend.text = element_text(size=16),
          legend.title = element_text(size=16))+
    xlab("Average Marginal Component Effects on Support for Policy Plan")+
    scale_color_grey(start = 0.1, end = 0.5)
  
  print(plgs)
  
  ggsave(filename = paste0("Graphs/weighted",filename ,".png"),
         plot = plgs,
         width = 12, height = 8)
  
  return(plgs)
  
}


################Figures for Paper:######

reg_all <- wtp_baseline_reg(data_wtp, "all")

plot1 <- plot_conjoint(reg_all, "Treatment Effects Policy Plan Conjoint, \n All Participants", "wtp")


#split sample plot

plot2 <- splitsampleplots_wtp(data_wtp, "highsupportintervention", "", "wtp_gs")+theme(axis.text.y = element_blank(),
                                                                                       legend.position = "right")


weighted_results <- ggpubr::ggarrange(plot1, plot2, nrow=1, legend="right", widths = c(0.5, 0.5))


ggsave(filename = paste0("Graphs/wtp_weighted_main_graphs.png"),
       plot = weighted_results,
       width = 20, height = 10)



###########Regression Tables for Paper########

#table 1
reg.wtp.1 <- wtp_baseline_reg(data_wtp, "all")

reg.wtp.2 <- wtp_baseline_reg_contr(data_wtp, "all")

reg.wtp.3 <- wtp_interactions(data_wtp, "all", "highsupportintervention")

reg.wtp.4 <- wtp_interactions_contr(data_wtp, "all", "highsupportintervention")

reg.wtp.5 <- wtp_interactions(data_wtp, "all", "supportintervention")

reg.wtp.6 <- wtp_interactions_contr(data_wtp, "all", "supportintervention")

row1 <- t(c("Demographic Controls",  "Yes", "Yes", "Yes")) %>% as.data.frame()


table.1 <- modelsummary(list(reg.wtp.2, reg.wtp.4, reg.wtp.6),
                        gof_omit="AIC|BIC", stars=T,
                        coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                        "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                        "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                        "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                        "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                        "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                        "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                        "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                        "competitionCarbon Border Tax" = "Carbon Border Tax",
                                        "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                        "highsupportintervention" = "Belief in St. Int. Bin.",
                                        "supportintervention"  = "Belief in St. Int. Cont."),
                        coef_omit = "region|age|gender|income", add_rows = row1,output = "latex")


header.t1 <- "DV: Support for Policy Plan"

tableheader.t1 <- c(c(" " = 1,  header.t1 = 3 ))

# set vector names 
names(tableheader.t1) <- c(" ", header.t1)



regtable.t1 <- table.1 %>% kableExtra::add_header_above(header=tableheader.t1) 


save_kable(regtable.t1, "Tables/wtp_weighted_conjoint_base.tex")



#table 2

reg.wtp.2.1 <- wtp_interactions_contr(data_wtp, "all", "highsupportinterventionecon")

reg.wtp.2.2 <- wtp_interactions_contr(data_wtp, "all", "supportinterventionecon")

reg.wtp.2.3 <- wtp_interactions_contr(data_wtp, "all", "highsupportinterventionfrail")

reg.wtp.2.4 <- wtp_interactions_contr(data_wtp, "all", "supportinterventionfrail")

reg.wtp.2.5 <- wtp_interactions_contr(data_wtp, "all", "highsupportinterventionworkers")

reg.wtp.2.6 <- wtp_interactions_contr(data_wtp, "all", "supportinterventionworkers")

table.2 <- modelsummary(list(reg.wtp.2.1, reg.wtp.2.3, reg.wtp.2.5),
                        gof_omit="AIC|BIC",
                        stars=T,
                        coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                        "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                        "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                        "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                        "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                        "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                        "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                        "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                        "competitionCarbon Border Tax" = "Carbon Border Tax",
                                        "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                        "highsupportinterventionworkers" = "Belief in St. Int. Bin. Workers",
                                        "supportinterventionworkers"  = "Belief in St. Int. Cont. Workers",
                                        "highsupportinterventionfrail" = "Belief in St. Int. Bin. Frail",
                                        "supportinterventionfrail"  = "Belief in St. Int. Cont. Frail",
                                        "highsupportinterventionecon" = "Belief in St. Int. Bin. Econ.",
                                        "supportinterventionecon"  = "Belief in St. Int. Cont. Econ."),
                        coef_omit = "region|age|gender|income", add_rows = row1,output = "latex")


workers <- "Just Workers"
frail <- "Just Sick and Elderly"
econ <- "Just Economy"

tableheader.tst2 <- c(" " = 1, econ=1, frail = 1, workers=1)
names(tableheader.tst2) <- c(" ", econ, frail, workers)

high <- "High"
low <- "Low"



regtable.t2 <- table.2 %>% kableExtra::add_header_above(header=tableheader.tst2) %>% kableExtra::add_header_above(header=tableheader.t1) 


save_kable(regtable.t2, "Tables/wtp_weighted_conjoint_diff_gs.tex")

#plot different interactions split sample

int.1  <- splitsampleplots_wtp(data_wtp, "highsupportinterventionecon", "\n Just Economy", "wtp_gs_w") + xlab("")

int.2  <- splitsampleplots_wtp(data_wtp, "highsupportinterventionfrail", "\n Just Sick and Elderly", "wtp_gs_w")+theme(axis.text.y = element_blank())+ xlab("AMCE on Support for Policy Plan")

int.3  <- splitsampleplots_wtp(data_wtp, "highsupportinterventionworkers", "\n Just Workers", "wtp_gs_w")+theme(axis.text.y = element_blank()) + xlab("")

intwtp.all <- ggpubr::ggarrange(int.1, int.2, int.3, nrow=1, common.legend = TRUE, legend="bottom", widths = c(0.8, 0.5,0.5))
intwtp.all

ggsave(filename = "Graphs/wtp_weighted_conjoint.indintsplit.png",
       plot = intwtp.all,
       width = 22, height = 10)

#table 3: only those who paid attention

reg.wtp.1v <- wtp_baseline_reg_contr(data_wtp, "wtpattention")

reg.wtp.2v <- wtp_interactions_contr(data_wtp, "wtpattention", "highsupportintervention")

reg.wtp.3v <- wtp_interactions_contr(data_wtp, "wtpattention", "supportintervention")


table.1a <- modelsummary(list(reg.wtp.1v, reg.wtp.2v, reg.wtp.3v),
                         gof_omit="AIC|BIC",
                         stars=T,
                         coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                         "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                         "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                         "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                         "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                         "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                         "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                         "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                         "competitionCarbon Border Tax" = "Carbon Border Tax",
                                         "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                         "highsupportintervention" = "Belief in St. Int. Bin.",
                                         "supportintervention"  = "Belief in St. Int. Cont."),
                         coef_omit = "region|age|gender|income",
                         output = "latex", add_rows = row1)


tableheader.t2 <- c(c(" " = 1,  header.t1 = 3 ))

# set vector names 
names(tableheader.t2) <- c(" ", header.t1)

regtable.t1a <- table.1a %>% kableExtra::add_header_above(header=tableheader.t2) 


save_kable(regtable.t1a, "Tables/wtp_weighted_conjoint_attention.tex")

reg_att <- wtp_baseline_reg(data_wtp, "wtpattention")

plota <- splitsampleplots_wtp(data_wtp[data_wtp$wtpattention==1,], "highsupportintervention", "\n Passed Attention Check Only", "wtp_attention")



#table 4: robustness check to other preferences

reg.wtp.1r <- wtp_interactions_robust(data_wtp, "all",  "highsupportintervention")

reg.wtp.2r <- wtp_interactions_robust(data_wtp, "all",  "supportintervention")

rowr1 <- t(c("Interaction Controls Other Attitudes", "Yes" , "Yes")) %>% as.data.frame()


table.r <- modelsummary(list(reg.wtp.1r , reg.wtp.2r),
                        gof_omit="AIC|BIC", stars=T,
                        coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                        "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                        "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                        "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                        "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                        "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                        "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                        "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                        "competitionCarbon Border Tax" = "Carbon Border Tax",
                                        "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                        "highsupportintervention" = "Belief in St. Int. Bin.",
                                        "supportintervention"  = "Belief in St. Int. Cont."),
                        coef_omit = "party|climate_policyn|climatehumans|left_right|region|age|gender|income", add_rows=rowr1
)



tableheader.tr <- c(c(" " = 1,  header.t1 = 2 ))

# set vector names 
names(tableheader.tr) <- c(" ", header.t1)

regtable.r <- table.r %>% kableExtra::add_header_above(header=tableheader.tr) 


save_kable(regtable.r, "Tables/wtp_weighted_conjoint_robust.tex")

table.rf1 <- modelsummary(list(reg.wtp.1r), stars=T,
                          gof_omit="AIC|BIC",
                          coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                          "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                          "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                          "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                          "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                          "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                          "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                          "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                          "competitionCarbon Border Tax" = "Carbon Border Tax",
                                          "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                          "highsupportintervention" = "Belief in St. Int. Bin.",
                                          "supportintervention"  = "Belief in St. Int. Cont.",
                                          "climatehumans" = "Belief in MM CC",
                                          "Climate_policyn" = "Support Climate Pol.",
                                          "party_2021_2" = "2021 Party Vote"),
                          coef_omit = "region|age|gender|income",
                          estimate  = "{estimate} ({std.error})",
                          statistic = NULL,
                          output="latex")

tableheader.tr1 <- c(c(" " = 1,  header.t1 = 1 ))

# set vector names 
names(tableheader.tr1) <- c(" ", header.t1)

regtable.rf1 <- table.rf1 %>% kableExtra::add_header_above(header=tableheader.tr1) 


save_kable(regtable.rf1, "Tables/wtp_weighted_conjoint_robust_full1.tex")


table.rf2 <- modelsummary(list(reg.wtp.2r), stars=T,
                          gof_omit="AIC|BIC",
                          coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                          "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                          "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                          "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                          "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                          "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                          "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                          "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                          "competitionCarbon Border Tax" = "Carbon Border Tax",
                                          "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                          "highsupportintervention" = "Belief in St. Int. Bin.",
                                          "supportintervention"  = "Belief in St. Int. Cont.",
                                          "climatehumans" = "Belief in MM CC",
                                          "Climate_policyn" = "Support Climate Pol.",
                                          "party_2021_2" = "2021 Party Vote"),
                          coef_omit = "region|age|gender|income",
                          estimate  = "{estimate} ({std.error})",
                          statistic = NULL,
                          output="latex")



regtable.rf2 <- table.rf2 %>% kableExtra::add_header_above(header=tableheader.tr1) 


save_kable(regtable.rf2, "Tables/wtp_weighted_conjoint_robust_full2.tex")




#################marginal means####################
data_wtp$competition <- as.character(paste(data_wtp$competition))
data_wtp$competition[data_wtp$competition=="No Change"] <- "No Change Comp"
data_wtp$competition <- as.factor(paste(data_wtp$competition))

mms <- cj(data_wtp, choice_plan ~ cost + effectiveness + compensation + competition, id = ~ResponseId, 
          estimate="mm")

mmstable <- kable(mms)

save_kable(mmstable, "Tables/wtp_weighted_mm.tex")


plot.data.mms <- data.prep.wtp.mm(mms)

pl.mm <- ggplot(data=plot.data.mms)+
  geom_point(aes(x = estimate, y = factor))+
  geom_segment(aes(x=lower, xend = upper,
                   y = factor, yend = factor))+
  geom_vline(xintercept = 0.5, lty = "dotted")+
  theme_bw()+
  scale_y_discrete(NULL, drop=FALSE,
                   labels = rev(c("Effectiveness",
                                  "25% CO2 Reduction by 2030",
                                  "55% CO2 Reduction by 2030",
                                  "Cost",
                                  "150EUR per Year",
                                  "300EUR per Year",
                                  "600EUR per Year",
                                  "1200EUR per Year",
                                  "Compensation Policy",
                                  "No Change",
                                  "Low Lump Sum Carbon Dividend",
                                  "High Lump Sum Carbon Dividend",
                                  "Low Progressive Carbon Dividend",
                                  "High Progressive Carbon Dividend",
                                  "Competitiveness Policy",
                                  "No Change",
                                  "Carbon Border Tax",
                                  "Environmental Clauses in PTAs")))+
  theme(axis.text.y = element_text(
    size = 16,
    face = rev(c('bold', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain','plain',
                 'bold', 'plain', 'plain', 'plain'))
  ),
  text = element_text(size=16),
  plot.title = element_text(size=18, hjust = 0.5),
  axis.title.x=element_text(size=16)
  )+
  xlab("Marginal Means of Support for Policy Plan")+
  ggtitle("Marginal Means Policy Plan Conjoint,
           All Participants")




########marginal means subgroups##########

data_wtp1 <- data_wtp %>% filter(highsupportintervention==1)
data_wtp0 <- data_wtp %>% filter(highsupportintervention==0)

mms1 <- cj(data_wtp1, choice_plan ~ cost + effectiveness + compensation + competition, id = ~ResponseId, 
           estimate="mm")

mms0 <- cj(data_wtp0, choice_plan ~ cost + effectiveness + compensation + competition, id = ~ResponseId, 
           estimate="mm")


mmstable_het_high <- kable(mms1)

save_kable(mmstable_het_high, "Tables/wtp_weighted_mm_het_high.tex")


mmstable_het_low <- kable(mms1)

save_kable(mmstable_het_low, "Tables/wtp_weighted_mm_het_low.tex")


plot.data.mms1 <- data.prep.wtp.mm(mms1)
plot.data.mms0 <- data.prep.wtp.mm(mms0)


pl.mm.het <- ggplot(data=plot.data.mms1)+
  geom_point(aes(x = estimate, y = factor,  color = paste("High", " n = ", nrow(data_wtp1)/10, sep=""), shape = paste("High", " n = ", nrow(data_wtp1)/10, sep="")))+
  geom_segment(aes(x=lower, xend = upper,
                   y = factor, yend = factor, color = paste("High", " n = ", nrow(data_wtp1)/10, sep="")))+
  geom_point(data= plot.data.mms0, aes(x = estimate, y = factor,  color = paste("Low", " n = ", nrow(data_wtp0)/10, sep=""), shape = paste("Low", " n = ", nrow(data_wtp0)/10, sep="")), position = position_nudge(y = -0.2))+
  geom_segment(data= plot.data.mms0, aes(x=lower, xend = upper,
                                         y = factor, yend = factor, color = paste("Low", " n = ", nrow(data_wtp0)/10, sep="")), position = position_nudge(y = -0.2))+
  geom_vline(xintercept = 0.5, lty = "dotted")+
  theme_bw()+
  scale_y_discrete(NULL, drop=FALSE,
                   labels = rev(c("Effectiveness",
                                  "25% CO2 Reduction by 2030",
                                  "55% CO2 Reduction by 2030",
                                  "Cost",
                                  "150EUR per Year",
                                  "300EUR per Year",
                                  "600EUR per Year",
                                  "1200EUR per Year",
                                  "Compensation Policy",
                                  "No Change",
                                  "Low Lump Sum Carbon Dividend",
                                  "High Lump Sum Carbon Dividend",
                                  "Low Progressive Carbon Dividend",
                                  "High Progressive Carbon Dividend",
                                  "Competitiveness Policy",
                                  "No Change",
                                  "Carbon Border Tax",
                                  "Environmental Clauses in PTAs")))+ 
  theme(axis.text.y = element_text(
    size = 16,
    face = rev(c('bold', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain','plain',
                 'bold', 'plain', 'plain', 'plain'))
  ),
  text = element_text(size=16),
  plot.title = element_text(size=18, hjust = 0.5),
  legend.position = "bottom",
  legend.title = element_text(size=16), #change legend title font size
  legend.text = element_text(size=16),
  axis.title.x=element_text(size=16)
  )+
  xlab("Marginal Means of Support for Policy Plan")+
  ggtitle("Marginal Means Policy Plan Conjoint,
           by Belief in State Intervention")+
  guides(color=guide_legend(title="Belief in State Intervention"),
         shape = guide_legend(title="Belief in State Intervention"))+
  scale_color_grey(start = 0.1, end = 0.5)




#########difference in marginal means###

data_wtp$highsupportintervention <- as.factor(data_wtp$highsupportintervention)

mm_diff <- cj(data_wtp, choice_plan ~ cost + effectiveness + compensation + competition, id = ~ResponseId, 
              estimate="mm_diff",  by = ~highsupportintervention)


mmstable_het_diff <- kable(mm_diff)

save_kable(mmstable_het_diff , "Tables/wtp_weighted_mm_het_diff.tex")

plot.data.mmdiff <- data.prep.wtp.mm(mm_diff)


pl.mm_diff <- ggplot(data=plot.data.mmdiff)+
  geom_point(aes(x = estimate, y = factor))+
  geom_segment(aes(x=lower, xend = upper,
                   y = factor, yend = factor))+
  geom_vline(xintercept = 0, lty = "dotted")+
  theme_bw()+
  scale_y_discrete(NULL, drop=FALSE,
                   labels = rev(c("Effectiveness",
                                  "25% CO2 Reduction by 2030",
                                  "55% CO2 Reduction by 2030",
                                  "Cost",
                                  "150EUR per Year",
                                  "300EUR per Year",
                                  "600EUR per Year",
                                  "1200EUR per Year",
                                  "Compensation Policy",
                                  "No Change",
                                  "Low Lump Sum Carbon Dividend",
                                  "High Lump Sum Carbon Dividend",
                                  "Low Progressive Carbon Dividend",
                                  "High Progressive Carbon Dividend",
                                  "Competitiveness Policy",
                                  "No Change",
                                  "Carbon Border Tax",
                                  "Environmental Clauses in PTAs")))+ 
  theme(axis.text.y = element_text(
    size = 16,
    face = rev(c('bold', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain',
                 'bold', 'plain', 'plain', 'plain', 'plain','plain',
                 'bold', 'plain', 'plain', 'plain'))
  ),
  text = element_text(size=16),
  plot.title = element_text(size=18, hjust = 0.5),
  axis.title.x=element_text(size=16)
  )+
  xlab("Difference in Marginal Means of Support for Policy Plan")+
  ggtitle("Difference in Marginal Means Policy Plan Conjoint,
           All Participants")




###merged figures for paper appendix
###marginal means:


pl.mm.het.2<- pl.mm.het  + theme(axis.text.y=element_blank(),
                                 axis.title.x= element_text(size=12),
                                 plot.title = element_text(size=14, hjust = 0.5))

pl.mm_diff.2 <- pl.mm_diff + theme(axis.text.y=element_blank(),
                                   axis.title.x= element_text(size=12),
                                   plot.title = element_text(size=14, hjust = 0.5))


pl.mm.2 <- pl.mm + theme(
  axis.title.x= element_text(size=12),
  plot.title = element_text(size=14, hjust = 0.5))

mmggplotwtp.all <- ggpubr::ggarrange(pl.mm.2, pl.mm.het.2, pl.mm_diff.2, nrow=1, common.legend = TRUE, legend="bottom", widths = c(0.8, 0.5,0.5))

ggsave(filename = "Graphs/wtp_weighted.mmall.png",
       plot = mmggplotwtp.all,
       width = 22, height = 10)


###table for split sample analyses#############

split_high <- wtp_baseline_reg_contr(data_wtp, "highsupportintervention")

split_low <- wtp_baseline_reg_contr(data_wtp, "lowsupportintervention")


split_high_a <- wtp_baseline_reg_contr(data_wtp[data_wtp$wtpattention==1,], "highsupportintervention")

split_low_a <- wtp_baseline_reg_contr(data_wtp[data_wtp$wtpattention==1,], "lowsupportintervention")

row2 <- t(c("Demographic Controls",  "Yes", "Yes", "Yes", "Yes")) %>% as.data.frame()


table_split <- modelsummary(list(split_high, split_low, split_high_a, split_low_a),
                            gof_omit="AIC|BIC",
                            stars=T,
                            coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                            "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                            "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                            "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                            "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                            "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                            "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                            "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                            "competitionCarbon Border Tax" = "Carbon Border Tax",
                                            "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs"
                            ),
                            coef_omit = "region|age|gender|income",
                            add_rows = row2,
                            output="latex")


tableheader.ts1 <- c(c(" " = 1,  header.t1 = 4 ))

# set vector names 
names(tableheader.ts1) <- c(" ", header.t1)

high_a <- "Passed Attention"
all <- "All"

tableheader.ts2 <- c(" " = 1, all=2, high_a = 2)
names(tableheader.ts2) <- c(" ", all, high_a)

high_b <- "High Belief St. Int."
low_b <- "Low Belief St. Int."

tableheader.ts3 <- c(" " = 1, high_b =1, low_b = 1, high_b=1, low_b=1)
names(tableheader.ts3) <- c(" ", high_b, low_b, high_b, low_b)



regtable.ts1 <- table_split %>% kableExtra::add_header_above(header=tableheader.ts3)%>% kableExtra::add_header_above(header=tableheader.ts2) %>% kableExtra::add_header_above(header=tableheader.ts1) 


save_kable(regtable.ts1, "Tables/wtp_weighted_conjoint_split_samples.tex")



####tables by treatment group
split_control <- wtp_baseline_reg_contr(data_wtp, "control")
split_compensationt <- wtp_baseline_reg_contr(data_wtp, "compensationt")
split_carbontax <- wtp_baseline_reg_contr(data_wtp, "carbontax")
interaction_treatments <- wtp_interactions_contr(data_wtp, "all", "vignettetreatment")

table_split_treatments<- modelsummary(list(split_control, split_compensationt, split_carbontax, interaction_treatments),
                                      stars=T,
                                      gof_omit="AIC|BIC",
                                      coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                                      "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                                      "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                                      "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                                      "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                                      "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                                      "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                                      "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                                      "competitionCarbon Border Tax" = "Carbon Border Tax",
                                                      "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                                      "vignettetreatment2" = "Compensation Households Vign.",
                                                      "vignettetreatment3" = "Carbon Tax Vign."
                                      ),
                                      coef_omit = "region|age|gender|income",
                                      add_rows = row2,
                                      output="latex")


tableheader.tvt <- c(c(" " = 1,  header.t1 = 4 ))

# set vector names 
names(tableheader.tvt) <- c(" ", header.t1)

control <- "Control"
comp <- "Compensation Households"
carbont <- "Carbon Tax"
interactions <- "Interactions"

tableheader.tvt2 <- c(" " = 1, control=1, comp = 1, carbont=1, interactions=1)
names(tableheader.tvt2) <- c("Vignette Treatment", control, comp, carbont, interactions)


regtable.tvt1 <- table_split_treatments%>% kableExtra::add_header_above(header=tableheader.tvt2) %>% kableExtra::add_header_above(header=tableheader.tvt) 

save_kable(regtable.tvt1, "Tables/wtp_weighted_conjoint_split_vignette_treatments.tex")




######plot of treatment groups split###
dg1 <- data.prep.wtp(split_control)
dg1$Treatment_Group <- "Control"

dg2 <- data.prep.wtp(split_compensationt)
dg2$Treatment_Group <- "Compensation"

dg3 <- data.prep.wtp(split_carbontax)
dg3$Treatment_Group <- "Carbon Tax"


plot.data.all.tg <- rbind(dg1, dg2, dg3)

pl <- ggplot(data=plot.data.all.tg)+
  geom_point(aes(x = Estimate, y = factor, color = Treatment_Group, shape = Treatment_Group), position = position_dodge(width=  .5))+
  geom_linerange(aes(x = Estimate, xmin=confintlower, xmax = confintupper, color = Treatment_Group, linetype = Treatment_Group,
                     y = factor, yend = factor), position = position_dodge(width= .5))+
  geom_vline(xintercept = 0, lty = "dotted")+
  theme_bw()+
  scale_y_discrete(NULL, drop=FALSE,
                   labels = rev(c("Effectiveness",
                                  "25% CO2 Reduction by 2030",
                                  "55% CO2 Reduction by 2030",
                                  "Cost",
                                  "150EUR per Year",
                                  "300EUR per Year",
                                  "600EUR per Year",
                                  "1200EUR per Year",
                                  "Compensation Policy",
                                  "No Change",
                                  "Low Lump Sum Carbon Dividend",
                                  "High Lump Sum Carbon Dividend",
                                  "Low Progressive Carbon Dividend",
                                  "High Progressive Carbon Dividend",
                                  "Competitiveness Policy",
                                  "No Change",
                                  "Carbon Border Tax",
                                  "Environmental Clauses in PTAs")))+
  xlab("AMCE")+ guides(color=guide_legend(title="Treatment Group"),
                       linetype=guide_legend(title="Treatment Group"),
                       shape = guide_legend(title="Treatment Group"))+
  scale_color_grey() 


pl <- pl + theme(axis.text.y = element_text(
  face = rev(c('bold', 'plain', 'plain',
               'bold', 'plain', 'plain', 'plain', 'plain',
               'bold', 'plain', 'plain', 'plain', 'plain','plain',
               'bold', 'plain', 'plain', 'plain'))
))

pl <- pl + ggtitle("Policy Plan Conjoint Results, \n by Vignette Treatment Group")

pl <- pl + theme(text = element_text(size=16),
                 axis.text.y = element_text(size=16),
                 plot.title = element_text(size=18, hjust = 0.5),
                 axis.title.x=element_text(size=16),
                 legend.position = "bottom"
)+xlab("Average Marginal Component Effects on Support for Policy Plan")


print(pl)


ggsave(filename = paste0("Graphs/wtp_weighted_vignette_groups.png"),
       plot = pl,
       width = 12, height = 8)


#################

###interactions by treatment group

control_interaction_treatments <- wtp_interactions_contr(data_wtp, "control", "highsupportintervention")
compensationt_interaction_treatments <- wtp_interactions_contr(data_wtp, "compensationt", "highsupportintervention")
carbontax_interaction_treatments <- wtp_interactions_contr(data_wtp, "carbontax", "highsupportintervention")



table_split_treatments_belief<- modelsummary(list(control_interaction_treatments, compensationt_interaction_treatments,carbontax_interaction_treatments),
                                             stars=T,
                                             gof_omit="AIC|BIC",
                                             coef_rename = c("cost1200 EUR per Year" = "Cost 1200 Eur per Year",
                                                             "cost600 EUR per Year" = "Cost 600 Eur per Year",
                                                             "cost300 EUR per Year" = "Cost 300 Eur per Year",
                                                             "effectiveness55% CO2 Reduction by 2030" = "55% CO2 Reduction",
                                                             "compensationHigh Lump Sum Carbon Dividend" = "High Lump Sum Carb. Div.",
                                                             "compensationHigh Progressive Carbon Dividend" = "High Progressive Carb. Div.",
                                                             "compensationLow Lump Sum Carbon Dividend" = "Low Lump Sum Carb. Div.",
                                                             "compensationLow Progressive Carbon Dividend"= "Low Progressive Carb. Div.",
                                                             "competitionCarbon Border Tax" = "Carbon Border Tax",
                                                             "competitionEnvironmental Clauses in PTAs" =  "Environmental Clauses in PTAs",
                                                             "highsupportintervention" = "Belief in St. Int. Bin."
                                             ),
                                             output="latex",
                                             coef_omit = "region|age|gender|income",
                                             add_rows = row1)



tableheader.tvt1 <- c(c(" " = 1,  header.t1 = 3 ))

# set vector names 
names(tableheader.tvt1) <- c(" ", header.t1)

tableheader.tvt2 <- c(" " = 1, control=1, comp = 1, carbont=1)
names(tableheader.tvt2) <- c("Vignette Treatment", control, comp, carbont)

regtable.tvt2 <- table_split_treatments_belief %>%  kableExtra::add_header_above(header=tableheader.tvt2) %>% kableExtra::add_header_above(header=tableheader.tvt1) 

save_kable(regtable.tvt2, "Tables/wtp_weighted_conjoint_split_vignette_treatments_beliefs.tex")



treat.1  <- splitsampleplots_wtp(data_wtp[data_wtp$control==1,], "highsupportintervention", "\n Control Group", "wtp_control") + xlab("")

treat.2  <- splitsampleplots_wtp(data_wtp[data_wtp$compensationt==1,], "highsupportintervention", "\n Compensation Treatment Group", "wtp_compensation") +theme(axis.text.y = element_blank())+ xlab("AMCE on Support for Policy")

treat.3  <- splitsampleplots_wtp(data_wtp[data_wtp$carbontax==1,], "highsupportintervention", "\n Carbon Tax Treatment Group", "wtp_control") + xlab("")+theme(axis.text.y = element_blank())


treatwtp.all <- ggpubr::ggarrange(treat.1, treat.2, treat.3, nrow=1, common.legend = TRUE, legend="bottom", widths = c(0.8, 0.5,0.5))
treatwtp.all

ggsave(filename = "Graphs/wtp_weighted_conjoint.treatmentsplit.png",
       plot = treatwtp.all,
       width = 22, height = 10)

